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We develop the kinetic theory of dilute cohesive granular gases in which the attractive part is 
described by a square well potential. We derive the hydrodynamic equations from the kinetic theory 
with the microscopic expressions for the dissipation rate and the transport coefficients. We check 
the validity of our theory by performing the direct simulation Monte Carlo. 
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I. INTRODUCTION 

The hydrodynamic description of granular materials is 
useful to know the rheological properties of the granular 
flow. Since granular materials are recognized to behave 
as unusual solids, liquids and gases, granular materials 
have attracted much interest among physicists [l|. The 
most idealistic granular system is a dilute gas without 
any external forces such as gravity. To analyze such a 
simple system is important to understand complex be¬ 
havior of granular materials. If the kinetic energy or the 
granular temperature of a granular gas homogeneously 
decreases because of inelastic collisions between grains, 
the time evolution of the temperature obeys Half’s law 
0. However, this homogeneous cooling state cannot be 
maintained as time goes on, because clusters of dense re¬ 
gion appear Such inhomogeneity of granular gases 

can be understood by granular hydrodynamics [6l4lOl| in 
which the transport coefficients for the inelastic hard core 
system for the dilute case [iMl and the moderately 
dense case [il[i3 can be determined by the inelastic 
Boltzmann-Enskog equation [^, [13, [T^ [l^ . These theo¬ 
retical results exhibit good agreements with the numeri¬ 
cal simulations, at least, for nearly homogeneous moder¬ 
ately dense granular flows [^. It should be noted 
that we often use the direct simulation Monte Carlo 
(DSMC) to evaluate the transport coefficients instead 
of using the molecular dynamics simulation, which was 
orig inally introduced by Bird [2^ to study rarefied gas 
|23l- l2^ and later has been extended to dilute inelastic 
gases^ [13 and to dense inelastic gases [H, H^l . This is 
because we should keep the system almost uniform. 

The interaction between contacting granular particles 
usually consists of the repulsive force and the dissipative 
force proportional to the relative speed. For fine powders 
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and wet granular particles, however, cohesive force can¬ 
not be ignored. The origins of such cohesive force are, 
respectively, van der Waals force for fine powders and 
capillary force for wet granular particles (SOl - IS^ - Such 
cohesive forces can cause the liquid-gas phase transition 
[sljl - the variations of cluster formation of freely falling 
granular particles and the enhancement of the 

jamming transition [^, j^. Thus, the study of cohe¬ 
sive granular materials is important for both physics and 
industry to treat real granular materials. In our previ¬ 
ous paper, we have demonstrated the existence of various 
phases for hne powders in the presence of a plane shear, 
which cannot be observed in granular gases under the 
shear [13. We have also developed the dynamic van der 
Waals model in describing such a system (sl] and obtain 
qualitatively consistent results with those in Ref. (d^ . 
These results suggest that the ordinary kinetic theory 
for a hard core system cannot be applied to this system. 
Needless to say, the kinetic theory is important to give us 
the microscopic basis of the macroscopic phenomenology 
such as Ref. [S^l and the simulation results such as Ref. 
[ 13 . In this paper, let us consider a granular gas whose 
interaction consists of the hard core for repulsive part 
and a square well potential for an attractive part. There 
exist some studies on the kinetic theory of gas molecules 
having the square well potential in which, the 

collision processes are categorized into four processes: (i) 
hard core collisions, (ii) entering processes, (iii) leaving 
processes from the well, and (iv) trapping processes by 
the well [d^, [d^, [d3 . Note that most of previous works 
study gases without dissipations in collisions except for 
some recent papers [13 [5l|> which do not discuss the 
transport coefficients. It should also be noted that some 
papers developed the kinetic theory based on different 
models for cohesion [H, [13 • 

In this paper, we derive modified Haff’s law and derive 
the transport coefficients for the dilute cohesive granular 
gases in freely cooling processes. For this purpose, we ex¬ 
tend the kinetic theory for the inelastic hard core system 
to the nearly elastic granular gases having the square well 
potential. The organization of this paper is as follows. In 
the next section, we evaluate the scattering angle for a 
two-body collision process as a function of the impact 
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parameter and the relative velocity of the colliding pair 
of particles by solving the Newton equation. In Sec. Hill 
we extend the kinetic theory for hard core granular gases 
to the gases having the square well potential to derive 
the transport coefficients in a set of the hydrodynamic 
equations. In Sec. lYl we compare them with those ob¬ 
tained by the DSMC. In Secs. lyland lVll we discuss and 
summarize our results, respectively. In Appendix we 
explain collision geometries for core collisions and grazing 
collisions to determine the velocity change during colli¬ 
sions in details. In Appendix m we briefly explain the 
procedure to obtain the transport coefficients by using 
the Chapman-Enskog theory. In Appendices [C] and |D1 
we calculate the second moment of the collision integral 
and two Sonine coefficients in terms of the kinetic theory, 
respectively. In Appendix |E1 we calculate the explicit ex¬ 
pressions of the transport coefficients in the high and low 
temperature limit. In Appendix |F1 we briefly summarize 
the DSMC algorithm. In Appendix [Gl we estimate the 
critical temperature, at which we cannot ignore the trap¬ 
ping process. 


II. SCATTERING ANGLE FOR THE SQUARE 
WELL POTENTIAL 

Let us calculate the scattering angle for monodisperse 
smooth inelastic hard spheres havin g th e square well po¬ 
tential whose mass is m [H, [53457l| . Here, the hard 

core potential associated with the square well attractive 
part for the relative distance r between two spheres is 
given by 


{ oo {r < d) 

—e {d < r < Xd) , (1) 

0 (r > Xd) 

where e and A are, respectively, the well depth and the 
well width ratio. We assume that collisions are inelastic 
only if particles hit the core (r = d) characterized by the 
restitution coefficient e. 





leave with the relative velocity v' after the scattering as 
depicted by Fig. [T] in the frame that the target is sta¬ 
tionary. The incident angle 6 between v and the normal 
unit vector k at the closest distance r = rmin between 
colliding particles is given by 


0 = b 


du 


0 ^l-b^u^-^U{l/u) 


( 2 ) 


where u= 1/r. Here, uq = l/rmin is the smaller one be¬ 
tween 1/d and the positive solution that the denominator 
of Eq. ([2]) is equal to zero [H, [s^, and k = is a 

unit vector parallel to ri 2 = r\ — r-i with the positions 
ri and r-i for particles 1 and 2, and ri 2 = |^’i 2 |- We have 
also introduced the impact parameter b for the incident 
process. Because the scattering is inelastic, in general, 
the impact parameter b' after the scattering and the an¬ 
gle 9' between k and v' differ from b and 9, respectively 
(Fig. [T]). Let us consider the case for b > Xd, where Eq. 
m reduces to 


nl/b 


9 = b 


du 


'o Vi — 2 


( 3 ) 


under the condition uq = 1/d. Because the particles do 
not collide, 9' = 9, the scattering angle y is given by 

X = TT -29 = 0, sin I = 0. (4) 

Next, we consider the case for 6 < Ad in which Eq. © 
can be rewritten as 


nl/Ad 


du 


9=b 


= arcsin 1^1+^ 


du 


0 y/l-b^u^ Ji/xd _ b2y^2 


4e 

rrup 


du 


li/\d \/v'^ - b'^u'^ ’ 


( 5 ) 


where we have introduced v as 


^ = \/l + —(6) 

V mv^ 

and uq = min (1/d, I'/b) with the introduction of a func¬ 
tion min(a;, y) to select the smaller one between x and y. 
We note that v is related to the refractive index . 

For b > vd, uq is given by uq = vjb and this collision is 
called a grazing collision (d^ . Is^ . Is^ . From Eq. m , we 
rewrite 9 as 

n ^ ■ f ■ f b \ 

9 = — + arcsin — arcsin —— . (7) 

2 V Ad / V vXd ) 


FIG. 1: A schematic view of a collision process. The dotted 
line represents the outer edge of the attractive potential. 


Let us consider a scattering process in which two parti¬ 
cles approach from far away with relative velocity v and 


Because the particle does not hit the core, 9' should be 
equal to 9. Then, the scattering angle y is given by 


= = TT — 29 = 2 arcsin 


X = X 


un - 1 — 

I izAd 


2 arcsin ( 

Xd 


( 8 ) 
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Equation ([H), thus, can be rewritten as 
X 

sin — = sin 
2 

Note that this collision does not exist for \ < v. 

For b < vd, uq is given by uq = 1/d, and then the 
particles hit the core of the potential. From Eq. m , we 


obtain 0: 



In this case, the collision is inelastic, and thus, 9' is not 
equal to 9. From the conservation of the angular momen¬ 
tum bv = b'v', 9' is given by 


arcsm 




V vXd 


— arcsm 


Ad 


(9) 


0' = arcsin 


= arcsin 



-f arcsin 
+ arcsin 



— arcsin 

— arcsin I 



e 




y/X^d? - 62 


+ 


b 


_^i 


cos^ 0 -I- O(e^), 

( 11 ) 


where we have introduced 0 as 

y/v'^d? — b'^ 


cos 0 = 


vd 


( 12 ) 


(see Appendix El for the derivation) and e = 1 — e. Thus, 
we obtain the scattering angle x as 


X 


= TT - 9 - 9'= + 0{e^) (13) 


with 


= tt — 2 arcsin ( -^i 
\XdJ 


xW = 


— 2 arcsin 
bv^ 


2 arcsin 
b 


V^^Adj 


.VA2d2 - 62 _ 52 

b 


y/XXXW^. 

We can rewrite Eq. (HI as 


• ^ ■ x^°^ , 1 ( 1 ) 


(14) 


cos"^ 0. (15) 


(16) 


TABLE I: Parameters corresponding to Fig. [2] 



(a) hard core 
(inelastic) 

(b) grazing 
(elastic) 

(c) no-collision 

6 

h/d < mm{v, A) 

min(:/, A) < b/d < A 

b/d> X 

. X 

sm 2 

Ea.(fT6l) 

Eq.® 

Eq.(|4]) 



(a) 




(c) 


These results are consistent with the previous study in 
the elastic limit (e —>■ 1) [i^. We regard the grazing col¬ 
lision as a combination of (ii) entering and (iii) leaving 
processes from the well [i^ . We ignore the trapping pro¬ 
cess by the attractive potential in the elastic limit (i. e. 
e —>■ 0) because colliding particles against hard cores have 
positive energies and the most of rebounding particles 
have still positive energies. In other words, if the trap¬ 
ping process is relevant, the inelastic Boltzmann equation 
is no longer valid. Thus, through the analysis of the in¬ 
elastic Boltzmann equation we will discuss whether it can 
be used even for weakly inelastic cohesive granular gases. 
We summarize the above results in Fig. [5] and Table HI 


FIG. 2: Schematic views of dynamic processes between two 
adjacent particles. There exist three types: (a) collisions via 
the hard core potential (inelastic), (b) grazing collisions (elas¬ 
tic), and (c) no-collisions. 


III. KINETIC THEORY AND 
HYDRODYNAMIC EQUATIONS 

If we consider a dilute and weakly inelastic homoge¬ 
neous granular gas, we may use the inelastic Boltzmann 












































equation 


where /(/, /) is the collision integral 
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=/(/,/), ( 17 ) 


J 


HfJ) ^ J J c^^0(min(A,i/) - b)\vi2 ■ k\ [xeCr{x,Vi2)fi'r',v'l,t)f{r,v!^,t) - a{x,vi2)fir,Vi,t)f{r,V2,t)] 

+ y y dke(b-min{X,iy))\vi2-k\[a(x,Vi2)fi'r',v'l,t)f{r,v'^,t)-a{x,vi2)fir,Vi,t)f{r,V2,t)]. (18) 

I- 


Here we have introduced the step function 0(a:) = 1 for 
a; > 0 and Q{x) = 0 otherwise. Here Vi 2 = |t’i 2 | with 
1^12 = Vi — V 2 with the velocity Vi {i = 1,2) for i-th par¬ 
ticle, cr(x,r)i 2 ) is the collision cross section between i-th 
and j-th particles, and 6 = 6/d is a dimensionless collision 
parameter. The factor Xe is related to the Jacobian of the 
transformation between pre-collisional velocities v ”, i?" 
and the velocities after collision Vi , V 2 [H, [H, l6ll | . 
The first and second terms on the right-hand-side of Eq. 
m correspond to inelastic and elastic collisions, respec¬ 
tively. For the sake of later discussion, we explicitly write 


the relationship 

between {vi,V 2 

) and 

{V1,V2) 


Vi = 

V'{ + ^Av, V2 

= u" 

1 ^ 

- Av, 

2 ’ 

(19) 

with 





At> = -2 ^1 

1 2 COS^ 0 \ 

cos^e)' 

y{2 • 

k)k + 0(e^) 

( 20 ) 

for inelastic hard core collisions 

and 




Av = - 2 ( 1 ; (2 

• k)k 


( 21 ) 


for elastic grazing collisions (see Appendix for the 
derivation). From Eq. ((2(111 , the explicit form of the factor 
Xe is given by 


X, = l+2e^»^^ + 0(t") (22) 

cos^ 0 

for inelastic hard core collisions. It should be noted that 
Eq. (|2^ is consistent with 1/e^ for inelastic hard core 
potential [n [ii, 113, im, because this can be expanded 
as 1/e^ = 1-I-2e J-C>(e^) in the nearly elastic limit and v 
and 0 reduce to v ^ 1 and 0 —>■ d, respectively, in the 
hard core limit from Eqs. (jH) and (HU). 


A. Homogeneous freely cooling 

In this subsection, let us determine the velocity dis¬ 
tribution function f{v,t) in freely cooling granular gases 
based on the Boltzmann equation (II3- First, we expand 


the distribution function in terms of Sonine polynomials 

[H [H 0,113, mi as 




1 -t- ^ ( cnS(, 
i.=i 


/ mV'^ \ 

V2^y 


(23) 


where V = |V| = |u — [/|is the local velocity 
fluctuation from the flow velocity U{r,t), /m(E) = 
n(TO/27rT)^/^ exp(—TOE^/2r) is the Maxwellian at the 
temperature T and the number density n, and S(,{x) = 
is the Sonine polynomial: 


Se\x) = Y, 


(-i)fcr(j+f + i) , 

r{j + k + i){i - k)ik\ 


(24) 


with the Gamma function r(a;). The time evolution of 
the granular temperature, obtained by the product of the 
Boltzmann equation with mv\/2 and integrating over Vi, 
is written as 

where we have introduced the cooling rate for the homo¬ 
geneous gas 


^( 0 ) = (26) 
3 V m 

Here, AI 2 is the second moment of the dimensionless col¬ 
lision integral 

M2 = -ydcic?/(/(°\/(°)), (27) 

where we have introduced the dimensionless velocity Ci = 
vi/vT{t) with the thermal velocity VT{t) = \j2T{t)/m, 
the dimensionless collision integral = 

and the dimensionless distribu¬ 
tion function /(^^(c) = {v^/n)f^^^{v,t). After some ma¬ 
nipulation of Eq. ([27]), AI 2 can be rewritten as [13, [63| 

M2=-^ J dci j dc2 J dk\ci2- k\a{x,ci2) 

X f°Hci)f^°\c2)A[cl+cl] (28) 
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with (t(x, C 12 ) = cr(x, V 12 )Id? and (t>{c) = exp(-c^), 

and Aip{ci) = ip{c'-) — '(/'(c^). It should be noted that the 
density keeps constant and the flow velocity is zero in the 
homogeneous state. 


B. Hydrodynamic equations 


In this subsection, let us derive the transport coeffi¬ 
cients which appear in a set of hydrodynamic equations. 
Multiplying the Boltzmann equation m by I, Vi and 
mvij2 and integrating over Vi, we obtain the hydrody¬ 
namic equations 


— + V • {nU) = 0, (29) 

^+ U-VU +—V ■ P = 0, (30) 

ot mn 

BT 2 

— + U-VT+ — {P:VU + V-q) + CT = 0, (31) 


where n{r,t) is the density field, U{r,t) is the flow veloc¬ 
ity, and T{r,t) is the granular temperature. The pressure 
tensor P, the heat flux q, and the cooling rate C are, re¬ 
spectively, defined as 


We, thus, rewrite the Boltzmann equation (El) as 




(/'' 


= / 


( 0 ) gfW 4 . 
(^/(o) + j/(i)+ ...), 1^/(0) + 5J(1) + ^ 


The equation at the zeroth order of Eq. 


to 




■ (39) 
is reduced 

(40) 


From Eqs (IM1)-(EID, the zeroth order hydrodynamic 
equations are, respectively, given by 


9(0) 

dt 


n = 0, 


9(0) 

dt 


U = Q, 


at 


T = (41) 


which are equivalent to those obtained in the previous 
subsection for the homogeneous cooling state. The ze¬ 
roth order of the pressure tensor and the heat flux are, 
respectively, given by 

4°^=nrdy, q(o)=0. (42) 

The first-order Boltzmann equation becomes 


III 

1 dvDij 

{V)f{r,v,t) + nTS^j, 

(32) 

q = 

1 dvS{V)f{r,v,t), 

(33) 

c = - 

m 1 
SnTj 

dvv^lifj), 

(34) 


where D,j{V) = m{V,Vj - V^S^j/3) and S{V) = 
{mV^l2 — 5T/2)V. We adopt the constitutive equations 
at the Navier-Stokes order 


P = -p[y,Uj+ V, [/. - -dy V • 17 , (35) 


(36) 


q = —k'VT — fi^n, 



(43) 


The corresponding first-order hydrodynamic equations 
are, respectively, given by 


—n = -V • (nU), 

5 ( 1 ) 1 

— U = -U-VU -V(nr), 

at mn 

f){i) 2 

—T=-U-VT--TV-U-C^^^T, (44) 


where p is the hydrostatic pressure, ij is the shear viscos¬ 
ity, K is the thermal conductivity, and p, is the coefficient 
proportional to the density gradient. Throughout this 
paper, we have assumed that the equation of the state 
p = nT is held because we are interested in the behavior 
in the dilute limit, though this assumption might not be 
true if the granular temperature is sufficiently low. 

To obtain the transport coefficients, we adopt the 
Chapman-Enskog method [H, HO, HH . Here, we expand 
the distribution function around Eq. 1)23^ as 

/ = /(o)+5/(i) + ... (37) 


by a small parameter S corresponding to the gradients of 
the helds. Similarly, the time derivative of the distribu¬ 
tion function is expanded as 


a _ a?'> 

m ~ 


(38) 


where the hrst-order dissipation rate is defined by 

(/<”'./“>). (45) 

We note that becomes zero because of the parity 
of the integral (HSl) [l^ . 1^ . 1^ . We assume that the 
distribution function depends on time and space only 
via its moments: the density n, the average velocity U 
and the temperature T as = f?^'>[v\n,U,T]. Then 
we can rewrite the first-order equation (1431) as 

= iV-U-V-Vn) + ^ (^TV -U-V- Vt) 

+ ^-((v.v)^-Tvp), 


(46) 
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where 


J(1) (/(0),/(i)) =_j(/(0),/(i))_/(/(i),/(0)) . (47) 

From the form of the first-order equation (H51) . the solu¬ 
tion of this equation is expected to have the form 

/(i) =A-V\ogT + B-Vlogn + C,jVjU^, (48) 


respectively, where 17® is given by 


17 ® = / dci / dc2 / dkd{x,ci2){ci2 ■ k)(j){ci)(j){c2) 


^ + ^aeSt{c\) 


f=i 


S(c 2 )-A[S(ci) + S(c 2 ) 

(55) 


where the explicit forms of the coefficients A, B, and Cij 
are given in Appendix [Bias Eqs. (IB19I) . (IB20|) . and (IB12I) . 
respectively. The pressure tensor and the heat flux can 
be written as 

P}p=-V V • , (49) 

= — kVT — ^Vn. (50) 

Substituting / = and Eq. into Eq. (15^ . 

we obtain the differential equation for the shear viscosity 
r] with respect to T as 


where 17® is given by 

17® = / dci / dc2 / dfca-(x,ci 2 )(ci 2 • fc)(/)(ci)0(c2) 


(51) 


^ + '^aeSe{cl] 


e=i 


Aj(c2)A Aj(ci) -b Aj(c2) 

(52) 


with Dij = Dij/e. Similarly, substituting Eq. o into 
Eq. (1331) , we obtain the differential equations for the ther¬ 
mal conductivity n and the coefficient /i with respect to 
T as 


A 

and 


(3C(°)Kr) + ^Kud^J^ni = (1 + 202) : 


2 m 


(53) 


- 3„Cl») A _ 3„,;(0) _ = A~, 


(54) 


with S = Syjm/e^. It should be noted that Eqs. (EU, 
(l5^ . and (l54l) are consistent with those in the previous 
study in the hard core limit [^ . 

C. Transport coefficients for the granular gases 
having the square well potential 


In the previous subsection, we have presented the gen¬ 
eral framework for the second moment (E51) and the dif¬ 
ferential equations of the transport coefficients dm), EH), 
and (IMl) in dilute granular cohesive granular gases with¬ 
out specification of mutual interactions between grains. 
In this subsection, let us derive the explicit forms of 
them for the square well potential outside and the hard 
core potential inside. Here, we assume that the zero¬ 
th order distribution function can be well reproduced by 
the truncation up to the third order Sonine polynomials 
[Hlellelil as 

/(°)(c) = (j){c) [l -b a 2 S 2 {c^) + a 3 S' 3 (c^)] , (56) 

where ai is automatically zero because the first order mo¬ 
ment is absorbed in the definition of the zeroth velocity 
distribution function. In this paper, we only consider the 
elastic limit e —^ 0. In addition, the coefficients 02 and 03 
can be, respectively, written as the series of e as shown 
in Appendix (Cl 


02 = 02 °^ -b ea' 2 '^ -b 0{e^) 
03 = 03 °^ -b + 0{e^) 


(57) 


where the coefficients are given by 


( 0 ) ( 0 ) 
ir, = a\ 


= 0 , 


HI) 


D ’ 


.( 1 ) 


D 


(58) 


with 


J 


poo POinax ^ ^ ^ / 1 \ P'^ __ __ 

N2 =2 J dci2 J — 6 ^)cf 2(5 — C12) exp 1 ^—-C12J J J b'— c'12) sin^ 


-12 


POO pOmax 

/ dci2 j 
/o Jo 


dbb{v^ - 6^)c!( 2(105 - 14c^2 “ cl2)exp ( -^c?2 


J dc [2 J db' 6'cf2(7 - c'll) sin^ exp , 


(59) 
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TVs =4 / dci 


db b{v - b )ci 2(105 - 14ci2 - 0 ^ 2 ) exp ( “ 2^12 


dc'i 2 / dfe' 6 'ci 2 sin^ exp ( —-c' 


dci 2 / - 6 ^) 012 ( 5 - 072 ) exp --C 12 


do'i 2 / db' b'c'^ 2 (7 + c'i|) sin^ x^°^' exp ( --c'i|) , 

(60) 


D = 


dci 2 J d6 6ci2sin^x^‘’^exp ^-ioi2^ y dc '^2 J - o'i2) sin^ x^°^'exp ^-^o'll^ 

j dci2 J db 6 oC2(7 - 0^2) sin^ exp J dc'^2 J db' 6 '0^2(7 + ofa) sin^ x^°^' exp ^-^o'll^ . 


(61) 


For simplicity we have introduced the notation = 
X^°\b',c',2). To obtain these expressions, we have ignored 
the terms proportional to 02 , a|, and 0203 because we are 
interested in nearly elastic situations. Therefore, from 
Eq. (1^ . we obtain 

M 2 =M^2^ + + Oie^), (62) 

where 


A4f = 0, (63) 



(64) 


with 6 max = min(z 2 (ci 2 ), A). Substituting Eqs. (1^51) and 
(1621) into Eq. (1251) . we obtain the time evolution of the 
temperature as the solid line in Fig. |31 in which the 
number density, the restitution coefficient, the potential 


width ratio, and the initial temperature are, respectively, 
nd^ = 0.05, e = 0.99, A = l.bd, and T = lOe. When we 
start from the temperature much higher than the well- 
depth, the decreases of the temperature obeys Haff’s law 
for hard core systems in the initial stage 0 . As the tem¬ 
perature approaches the well-depth, the rate of temper¬ 
ature decrease is larger than Haff’s law. A similar result 
on the crossover from Haff’s law to a faster decrease of 
the temperature has already been reported by Ref. 

Next, let us calculate the transport coefficients. Simi¬ 
lar to the previous case, with the dropping the contribu¬ 
tions from a|, a|, and 0203 , the coefficients and fl® 
defined in Eqs. (15^ and (|55]) are, respectively, given by 
(see Appendix iDl for the derivation) 

with 


dci2 y d6 60^2 sin^x^°^ exp , 

y dci 2 J (6660^2 ( 63 -18c?2 + Ci2)sin2 x'”^ exp 

- 4 ^^ ^^12 (693 - 2970^2 + 33c^ 2 “ ^2) sin^ X^°^ exp ic^2^ 

- y dci2 J db bcl2X^^'* sin 2x^°^ exp (^“^'^12^ 

-h J dci2 J dbb{v‘^ - b^)cl2 f - - sin^ exp ^ 

^e(0) ^ J gj^2 ^(0) g^p ^-ici2^ , 


( 66 ) 


(67) 

( 68 ) 









J ^^12 J ~^^l 2 (63 - 18 c ?2 + 0 ^ 2 ) sin^ exp (^-^<^ 12 ^ 

+ 4 ^^ dci 2 J db 6 ci 2 (693 - 297ci2 + 33ci2 - C12) sin^ exp ^- 


"12 


— J dci2 J d6 6c^2X^^^sin2x^°^exp f--C 12 


■ 


j-oo pb„ 

/ dci 2 / 

/o Jo 
/*00 


^ ~ ~ / 1 
db b(v^ - &^)ci2 cos^ ^ exp ( - -C12 


\/^ 7°° /-Omax _ _ _ / 1 , 

H-— / (ici2 / c?& 5(1^^ — &^)cf2 (25 — IIC12) exp I —-c 

8 ./n ./n V 2 


( 69 ) 


It should be noted that the zeroth order of these quanti¬ 
ties, Eqs. (l66l) and (1681) . are the exactly same as the ones 
obtained by the previous study [43 |. 

Let us perturbatively solve the differential equation of 
the shear viscosity (ICT) with respect to the small param¬ 
eter e. We expand the shear viscosity as 

r; = 77(°)-bO(e2). (70) 

From Eqs. (ED), (EH), and dZOl), we rewrite the differential 
equation of the shear viscosity ED as 


— -ndrx — 


-\n(f 

5 

= nT. 


m 

m 




T 


±( 

dT \ 


,(o) 


£77 


( 1 ) 


— ( + e n<°'> +■■■ 


,(o) 


£77 


( 1 ) 


(71) 


Solving the zeroth and first order of this equation, we 
obtain 


77 ( 0 ) = 


77 ( 1 ) = 



(72) 

(73) 


Similarly, the thermal conductivity k and the coefficient 
/i are, respectively, given by 


K = K^^'> + + 0{e^), (74) 

= 7i(o) -H £7,(1) + c>(e2) (75) 


with 


K' ' = — 


(0) 

75 


1 


16d2 

V m 

(1) 


■k(°) 

75 


1 

0 


8d2 


5 

1 



2T a 




( 1 ) 


m ^ 


e(0) 


dT 




d 


(°) = 0, 


^(1) _ 5 7W^^^k(°)T 


2n 




(0) 



8 nd'^ V 2m J7®(°) 


(76) 


(77) 

(78) 

(79) 


We note that the zeroth order terms of these transport 
coefficients, Eqs. dZH) and dZH) are identical to those ob¬ 
tained by the previous studies (i^ . 

We obtain the expressions of the transport coefficients 
as Eqs. ED, dZQl), ED, and ED. The above procedure is 
not practically useful for the simulation of the hydrody¬ 
namic equations because we need to calculate the double 
integrals at every step. To reduce the calculation cost, 
we compare the results with high and low temperature 
expansions. From the calculation in Appendix |E1 we can 
obtain the explicit expressions of the dissipation rate and 
the transport coefficients as in Table|lTl As a final remark 
in this section, we note that our results up to 02 order 
in Eq. (1561) are almost identical to those up to 03 in the 
elastic limit. This ensures that the expansion around the 
Maxwellian gives well converged results by Eq. (IMl) . 


J 


IV. COMPARISON WITH THE NUMERICAL 
RESULTS 

To check the validity of the kinetic theory, we com¬ 
pare the transport coefficients derived from the kinetic 


theory in the previous section with those obtained by 
the DSMC, which is known as the accurate numerical 
method to solve the Boltzmann equation |65| . We 

note that stochastic treatment of collisions via DSMC 
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TABLE II: High temperature expansion of each quantity and low temperature expansion of the second moment up to first 
order of e/T and e. 


Afz = 2v/^e (^1 + (T-^oo), M 2 = 2V^e (^1 + \^(T ^ 0) 

L!:; = -4v^[l + e 
fij = —4%/^ 

V 


1 + e 


16d2 

75 


mT 

TT 

T 


ti = e 


64:(P V vrm 
1185 


1280 T 96 

1989 e A - 1 

1280 ~ T 96 t 

1567 e A - 1 

l + eTTTTTTT + 


2(15A'‘ + 15A® + 2A^ + 2A + 2) + 3A^(A + 1)(5A^ - 1) log ^ 

A + 1 

A- 1 


'3840 T 96 { 


2(15A'‘ + 15A® + 2A^ + 2A + 2) + 3A^(A + 1)(5A^ - 1) log 

I 


2(15A'‘ + 15A® + 2A^ + 2A + 2) + 3A^(A + 1)(5A^ - 1) log 


A + 1 
A- 1 


, 539 £ A - 1 

l + eT7t7t7t + 


1280 T 96 { 


2(15A'‘ + 15A® + 2A^ + 2A + 2) + 3A^(A + 1)(5A^ - 1) log 


A + 1 
A- 1 
A + 1 


1024nd2 V vrm 



FIG. 3: (Color online) The time evolution of the granular 
temperature for ndP = 0.05, A = 1.5, and e = 0.99 obtained 
by the kinetic theory (blue solid line) and that by the DSMC 
(red open circles), where t* = t^/Tfrajd and the initial tem¬ 
perature is set to be lOe. The dotted line represents Haft’s 
law for inelastic hard core spheres in which each particle has 
the diameter d. 


ensures the system uniform, which is suitable to measure 
the transport coefficients. 


A. Cooling coefficient 

In this subsection, we check the time evolution of the 
granular temperature for homogeneous cooling state and 
the second moment Ad 2- We prepare monodisperse N 
particles in a cubic box with the linear system size L. 
We distribute particles at random as an initial condition, 
where the initial velocity distribution obeys Maxwellian 
with the temperature T = lOe. Figure [3] shows the 
time evolution of the temperature obtained by the DSMC 
and Eq. (1351) . in which the number of particles, the sys¬ 
tem size, the number density, the potential width, and 
the restitution coefficient are, respectively, N = 6, 250, 
L = 50d, nd^ = 0.05 A = 1.5, and e = 0.99. The time 
evolution obtained by the kinetic theory fairly agrees 


with that by the DSMC. Figure S] shows the compari¬ 
son of the second moment Ad 2 obtained by the kinetic 
theory with that by the DSMC, which is also consistent 
each other, where Ad2 at high temperature limit is iden¬ 
tical to that for the hard core system with the diameter 
d. 



Tie 


FIG. 4: (Golor online) The granular temperature dependence 
of the second moment Ad 2 obtained by the DSMC (red open 
circles) and that by the kinetic theory up to as order (blue 
solid line), where T* is the dimensionless temperature defined 
by T* = T/e. The dotted line represents Ad 2 for the hard 
core system with the diameter d. The dashed (dot-dashed) 
line represents Ads obtained from the high (low) temperature 
expansion. 


B. Shear viscosity 

Let us compare the result of the shear viscosity by the 
kinetic theory with that by the DSMC in this subsec¬ 
tion. The particles are distributed at random and the 
velocity distribution satisfies Maxwellian at the initial 
condition. Then, we apply the shear with the aid of the 
Lees-Edwards walls at y = ±L/2, whose z-component 
is ±Fwaii- In the initial stage, the energy injection from 
shear is not balanced with the energy dissipation. Then, 
as time goes on, the system reaches a nonequilibrium 
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FIG. 5: (Color online) A schematic view of our setup to mea¬ 
sure the shear viscosity. The walls at y = {y = —Ll2) 
move to positive (negative) 2 -direction. 


steady state. In this stage, we calculate the shear viscos¬ 
ity defined by 


= - lim (80) 

t—^OC 'J 

where 7 is a bulk shear rate defined by the gradient 
of the flow velocity Uz and Pxy can be measured by 
the DSMC. To suppress the boundary effects, we mea¬ 
sure 7 in the range —L/4 < y < that is, 7 = 

{Uz\y=L/i — Uz\y=-L/A)/{L/2). Although the Newtonian 
shear viscosity should be measured by a relaxation pro¬ 
cess from the initial perturbation for the homogeneous 
cooling system [ 23 , lH, [o^l , this method is hard to mea¬ 
sure the shear viscosity in the low temperature region. 
It is also noted that the Newtonian viscosity is known 
to be identical to the steady state shear viscosity in the 
elastic limit [h^, which is the reason why we adopt the 
above setup. Figure |6] shows the comparison of the shear 
viscosity obtained by the kinetic theory with that by 
the DSMC, in which the number of particles, the sys¬ 
tem size, the number density, the potential width, and 
the restitution coefficient are, respectively, N = 10,000, 
L = 3 , OOOd, nd? = 0.01 A = 2.5, and e = 0.99. Similar 
to the case of Ad 2 , the shear viscosity obtained by the 
DSMC is identical to that obtained from the kinetic the¬ 
ory for hard-core systems with a particle diameter, d, in 
the high temperature limit. We cannot measure the shear 
viscosity for T < 10“^£ because the system is heated up 
by the shear even if we start from a lower temperature. 
The first order solution of the kinetic theory with respect 
to e also deviates from the zeroth order solution below 
this temperature, which suggests that the hydrodynamic 
description is no longer valid in this regime. This may 
correspond to the limitation of the inelastic Boltzmann 
equation, where the trapping processes cannot be ignored 
even in the elastic limit. 


C. Thermal conductivity 

Next, we compare the thermal conductivity obtained 
by the kinetic theory with that by the DSMC. Although 



* 



FIG. 6: (Golor online) Granular temperature dependence of 
the shear viscosity obtained by the DSMC (red open circles), 
that by the elastic kinetic theory (black solid squares in the 
previous study and black dashed line), and that by the 
kinetic theory (blue solid line), where 77 * is the dimensionless 
shear viscosity defined as rj* = yd?/me. The dotted line 
represents the shear viscosity for the hard core system of the 
diameter d. The dot-dashed line represents the shear viscosity 
obtained from the high temperature expansion. 


X 



FIG. 7: (Golor online) A schematic view of our setup to mea¬ 
sure the thermal conductivity. The temperature of the left 
(right) side wall is kept at Tl (Tr). 


the heat flux contains the term proportional to the den¬ 
sity gradient, we ignore its contribution because the term 
disappears in the elastic limit e —>■ 1 as in Eq. dZHl). To ob¬ 
tain the thermal conductivity from the DSMC, we solve 
the heat equation under a confined geometry shown in 
Fig. [a where the temperature at the left (right) wall at 
y = -LI2 {y = L/2) keeps Tp (Tr) [IMIl. In the steady 
state, because hydrodynamic variables depend only on y, 
the heat equation (1^ is reduced to 


i-A,, = cr = -kAt. ( 81 ) 


Let us nondimensionalize the quantities using the mass 
m, the system size L, and the well depth e as 


Tl 

n=—, y = Ly*, T = eT*, (82) 

P=^P-. = ( 83 ) 
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Thus, we rewrite the heat equation as 


dy*'^ 


0 = -3720-1/3 


(84) 


with 9 — T*3/2 and 7 ^ = {1 j\/2)p*'^M.^! k,'*. By multi¬ 
plying d9/dy* in both sides of Eq. (I5T)) and integrating 
the equation from y* = 0 to y*, we obtain 


dO , 1 


(85) 


where C is given by C = 9'^ + with 0o = 0|j/*=o 

and 9'q — d9/dy*\yt=o. Here, we consider the system 
that the temperature at y = — L /2 is lower than that at 
y = L/2, in which the plus sign is selected in Eq. (1551) . 
Under this condition, the solution of Eq. (I55|) has the 
following form 


y = 


, 1/3 
(0 _ 

27 


—Q\/P^ — 02 -I- ,52 arctan 


0 




+ \//32 — 1 — ^2 arctan 




( 86 ) 


where (3 = {{9'^+ l}i/2 and 0 = (0/0o)i/3. 

To obtain k! from the DSMC, we numerically evaluate 
7 from the comparison of the temperature profile (1861) 
with that by the DSMC in the range —L/5 < y < L/10 
as in Fig. |S1 It should be noted that we omit the data 
near the walls to suppress the boundary effects. Using 
the estimated 7 and the simulation results 0o, 0g, and 
A 42 in the homogeneous freely cooling, we estimate k! in 
terms of the DSMC. Here, the number of particles, the 
system size, the number density, the potential width, and 
the restitution coefficient are, respectively, N = 10,000, 
L = 3, OOOd, nd^ = 0.01 A = 2.5, and e = 0.99. Figure 
|9] shows the results of the DSMC and the kinetic theory, 
which is similar to that for rj. The heat conductivity 
in the high temperature limit of DSMC is identical to 
that predicted by the kinetic theory for hard-core systems 
with a particle diameter d. We note that the profile of the 
temperature described by Eq. (|55)) cannot be achieved 
for T < IQ-^e. Moreover, the perturbative contribution 
becomes larger than the base value of the perturbation 
dzi for T < O.le as in the case of the viscosity. 


V. DISCUSSION 

In this paper, we have obtained the transport coef¬ 
ficients as functions of the granular temperature. The 
transport coefficients in high temperature limit are iden¬ 
tical to those for the hard core system with the diameter 
d. Let us consider this reason. As explained in Sec. m the 
collision is inelastic for b < min(z/(i. Ad) while it becomes 
an elastic grazing collision for min(i/(i. Ad) < b < Ad. The 
value oi v = y^l -I- AeKmv^) converges to 1 in the high 



y/L 


FIG. 8 : (Color online) The solution of the heat equation (blue 
solid line) and the temperature profile obtained by the DSMC 
(red open circles). We choose 7 to fit the DSMC result in the 
range —L/5 < y < L/10. 



* 



FIG. 9: (Color online) The temperature dependence of the 
thermal conductivity obtained by the DSMC (red open cir¬ 
cles), that by the elastic kinetic theory (black solid squares 
in the previous study and black dashed line), and that 
by the kinetic theory (blue solid line), where k* is the dimen¬ 
sionless thermal conductivity defined as k* = Kd^y^m/e. The 
dotted line represents the thermal conductivity for the hard 
core system of the diameter d. The dot-dashed line repre¬ 
sents the shear viscosity obtained from the high temperature 
expansion. 


temperature limit. On the other hand, grazing collisions 
only change the directions of colliding particles and the 
kinetic energy is kept unchanged. Therefore, the energy 
change by collisions in high temperature limit is identical 
to that for the hard core system of the diameter d. 

Below T ~ 10“ ^e, the first order solutions of the trans¬ 
port coefficients with respect to e deviate from the ze¬ 
roth order solutions. Moreover, the first order solutions 
diverge as T~^ in the low temperature limit. This is 
because v diverges as 


v = 



rp-ll2 


(87) 


in the low temperature limit. This indicates that our 
hydrodynamic description in terms of the perturbation 
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method is no longer valid for low temperature. 

Murphy and Subramaniam [5l| studied the homoge¬ 
neous cooling state for a system of particles having an 
inelastic hard core associated with van der Waals po¬ 
tential. They obtained that the time evolution of the 
granular temperature obeys Haff’s law in the initial stage 
and decreases faster as time goes on, then approaches to 
Haff’s law for e = 0. They considered that the parti¬ 
cles aggregate after the collision when two particles have 
small kinetic energy with compared to the potential well 
keeping the potential contribution after the coalescence. 
Although we do not consider the aggregation process, the 
time evolution of the granular temperature in Fig. [3] is 
similar to their result. 

Our theory becomes invalid for T < O.le as shown in 
Figs, m and ini Let us estimate this critical temperature of 
coalescence processes from a simple one dimensional col¬ 
lision model. As explained in Appendix [U] if the kinetic 
energy is less than the well depth, the particle cannot 
escape from the well and be trapped by another parti¬ 
cle. This critical velocity can be estimated as Utrap — 
{8(1 — e)e/TO}^/^, which leads to the corresponding crit¬ 
ical temperature as Ttrap = ~ 4(1 — e)e. Us¬ 

ing our choice of parameter (e = 0.99), this temperature 
becomes Ttrap = O.Ode, which qualitatively reproduces 
the lower bound of our theory as shown in Figs. [Bland IHl 
Even if we can ignore aggregations of colliding particles, 
the equation of state p = nT is no longer valid for low 
temperature regime. The replacement of the equation 
of state will be discussed elsewhere. It should be noted, 
however, that realistic situations might not be described 
by Smoluchowski’s rate equation as used in Refs. 
because the biggest cluster may absorb other particles 
[43| . We will study the effects of aggregation processes 
in the near future. 

Here, we have only focused on the dilute system. To 
discuss the behavior of a system with finite density is also 
our future work. 


VI. CONCLUSION 

In this paper, we have developed the kinetic theory 
for dilute cohesive granular gases having the square well 
potential to derive the hydrodynamic equations using 
the Champan-Enskog theory for the inelastic Boltzmann 
equation. We have obtained the second moment M .2 of 
the collision integral and the transport coefficients for 
this system. We have found that they are identical to 
those for hard core gases at high temperature and the 
hydrodynamic description is no longer valid at low tem¬ 
perature. We have also performed DSMC simulation to 
check the validity of the kinetic theory and found that 
all results of DSMC are consistent with those obtained 
by the kinetic theory. 
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Appendix A: Collision geometry for the square well 
potential 



FIG. 10: Collision geometry for a grazing collision. Two parti¬ 
cles approach from Oi and leave for O2. The solid and dotted 
circles represent the hard core (radius d) and the outer edge 
of the potential (radius Ad), respectively. 

In this appendix, let us explain the collision geometry 
scattered by the square well potential. First, we consider 
the case for a grazing collision as in Fig. (TUI in the frame 
that the target is stationary. Let us consider the process 
that two particles approach from far away with relative 
velocity v from Oi. When the incident particle enters the 
well at the point A, the relative velocity changes because 
of the conservation of the energy and the angular mo¬ 
mentum, whose speed inside the well is given by vv. At 
the point A, the relative velocity perpendicular to OA is 
conserved, that is, usina = vvs\nj3 is satisfied [H^. The 
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change of the velocity parallel to OA is given by 


From these equations, we can rewrite Eq. as 


vv cos j3 — V cos a = vv\ll -^ sm' a — v cos a 


1 




= \ y — sin^ a — cos a 


) (Al) 


which means that the velocity change Ava at the point 
A satisfies 


Aua = - (^\/i 


— sin^ a — cos a ) vva^ (A2) 


with the unit vector rA = (cos( 7 r — a), sin( 7 r — a))"*" par¬ 
allel to OA. We note that the minus sign in Eq. (IA2I) 
comes from the fact that the velocity change is opposite 
direction to Fa- 

Similarly, the component of the velocity change par¬ 
allel to OC at the point C is given by (cos a — 
\J— sin^ 0 !)v, which means that the velocity change 
Avc at the point C becomes 


Avc 



— sm a — cos a 


vrc 


(A3) 


A« = 2«cos«('“}'“S 

y sin( 7 r — 0) 

„ . /cos( 7 r - 0)\ 

= _2^;cOs(^-0) 

= —2{v ■ k)k, (AS) 

with the unit vector k = (cos( 7 r — 6 ), sin( 7 r — 0 ))"^. 



with the unit vector vq = (cos( 7 r -20-1- a), sin( 7 r -20-1- 
a))T. 

From Eqs. (IA2I) and (IA3I) . the velocity change Av dur¬ 
ing this grazing collision becomes 

Av = Ava + Avc 

= —2 (^\/v^ — sin^ a — cosa^ 

xucos(0-a)(^^;;^®y:^jy (A4) 


FIG. 11: Collision geometry for a core collision. Two particles 
approach from Oi and leave for 02.The solid and dotted lines 
represent the hard core (radius d) and the outer edge of the 
potential (radius Ad), respectively. 

Next, let us consider the case for a hard core collision 
as in Fig.[TT] In this case, an inelastic collision takes place 
at the point D. To calculate the energy dissipation at the 
point D, we consider the angle 0 between the relative 
velocity of the particle and OB. From AB = Adsin(0 —a), 
BD = OB — OD = (Acos(0 — a) — l)d, we can write 0 as 


From Eq. o and a = arcsin(AE/OA) = arcsin(6/Ad), 
the following relationships are satisfied: 


^ AD Asin(0 —a) 

tan0 =-=-- -A—. 

BD Acos(0 — a) — 1 


(A9) 


cos(^ — a) = cos ( — — arcsin —^ ^, (A5) 

\ 2 v\d) uXd 

a ' f ' ^ ^ 

cos a = sm arcsin —— — arcsin — 

\ vXd Xd 

. ( . b \ ( , h 

= sm arcsin —— cos arcsin — 

\ vXd J \ Xd 

( • M • ^ ^ 

— cos arcsin —— sm arcsin — 

\ vXd J \ Xd 

(\/A2d2 _ 62 _ _ 52^ ^ 

(A6) 


vX^d? 


and 


■\/^2 _ 


sm a — cos a 


^ (yv'^X^d? - 02 _ y'A2d2 _ (A7) 


From Eq. (nni, cos(0 —a) and sin(0 —a) are, respectively, 
given by 


cos 


(0 — a) = cos arcsin —- — arcsin —— 
ud vXd ^ 

v'^d? — b'^y/u'^X^d? — , 

(AlO) 


AAd2 


sin (0 


a) = sin 


( . h . b \ 

arcsin —- — arcsin —— 

\ vd vXdJ 


= {y/v'^X^d? — 02 — y/iX^ 

1^2 Ad2 V 


02), 

(All) 


and substituting Eqs. (jAlOl) and (IA11|1 into Eq. (|A9|1 . we 
obtain 


tan 0 = 


b 

y/v'^d'^ — 02 ’ 


(A12) 
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or, equivalently, Eq. (HU). From this, we can calculate 
the change Av'^ after the collision at the point B as 

Av'^ = —(1 — e^)iy^v^ cos^ 0 

= (A13) 

Correspondingly, the change of relative velocity Av is 
given by 


Av = — 

= -2 


(v ■ k) + y {v ■ fe)2 — (1 — cos^ 0 


k 


1 2 0 
2 ^ cos^ 0 


{v-k)k + 0{e‘^), (A14) 


which reduces to Av = —2(v ■ k)k in the elastic limit. 

Appendix B: Chapman-Enskog expansion 

In this Appendix, let us explain the outline of the 
Chapman-Enskog theory [l4l l60l|. As explained in Sec. 
uni the zeroth order distribution function, , is deter¬ 
mined by Eq. (HOI) in the form Eq. (1^ [^. The first 
order distribution satisfies Eq. (l46l) . which can be 
rewritten as 


,9(0) j(i) 


j(i) 


5/(0) 


dt y J ^ gj. 

= A-V\ogT + B-Vlogn + CijVjUi, (Bl) 

where the coefficients A, B, and Cij are, respectively, 
given by 

T d 




= V 


(y/(0)) 

T fmV^ 
m 


2T 


m dV 
VdV^ 


/(O) 


/(o), (B2) 


B(y) = -y/(o)_ 1^/(0) 
m oV 


(B3) 


C«(V) = A(i.;(o)) 


-Is,. 


= ( ^ 


d 

W 

1 




From Eq. (|Bip . is expected to have the form 

/(i) =A-V\ogT + B-V\ogn + CijWjU,. 

The relationships between the coefficients A, B, Cij and 
A, B, Cij are, respectively, obtained by substituting the 
solution Eq. (IB5I) into Eq. m as: 


1 

o 

+ 


=A, 

(B6) 



=B, 

(B7) 



=C^j, 

(B8) 


where we have used = 0 because the coefhcient Cij 
is traceless. 

Substituting Eq. (jB5l) into Eq. (15^ with the aid of Eqs. 
and dm, we obtain 


dVD,,iV)CkiiV)ViUk 


= -r]{w,U,+WjU, - -5,jV -U] . (B9) 


Therefore, the shear viscosity r] is given by 

ri=-^JdVD,,{V)C,,{V). (BIO) 

Substituting Eq. (I56l) into Eq. (IB4|) . we obtain the ex¬ 
plicit form of Cij{V) as 




-^D,,{V) |l + E I /m(C). 

(Bll) 


This form and Eq. (IB8I) leads to 


C.,(V) = ^A,(F)/m(E), 


(B12) 


where Ci is a constant. Substituting Eq. (IB12I) into Eq. 
(IBIOI) . we obtain Ci = —r]/{nT). 

Similarly, substituting f^^'> into Eq. (1551) with the aid 
of Eqs. and (1501) . we obtain 



[T j 

(B4) 



\n. 

(B5) 

Therefore, 


tivity and the coefficient ji as 


n = -— dVSiV)-AiV), (B15) 


M = / dVS{V)-B{V). 


(B16) 


Substituting Eq. (1501) into Eqs. (|B2I) and (IB3|) . we obtain 
the explicit forms of A{V) and B{'V) as 
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A(y) =y [l + ^ 

B{V) =Y,cieVSf^Jp{c^)fM{V). 


+ J2ae s[^^^\c^)Se{c^) + {l-c^)Sl^_!^\c^) \ fuiV), (B17) 


e =3 


(B18) 


Equations (IB6I) and (lB7l) leads to 

A = -^SiV)fM{V), (B19) 

S = -yW/m(F), (B20) 

where Ai and Bi are constants. Substituting Eqs. 112]) 
and (jB3p into Eq. (jB15p and (|B16I) . respectively, and 
integrating over V, we obtain Ai = 2mK/5nT and Bi = 
2mfil5T^. 

Let us determine the explicit forms of the transport co¬ 
efficients. Multiplying Eq. (IB8I1 by Dij{Vi) and integrate 
over Vi, we obtain 

ioc(o)t^ + J dViA,(Vi)jW 

= J dViD,j{Vi)a,{Vi). (B21) 

The second term on the left-hand-side of Eq. (IB21I) is 
written as 

j (B22) 

where 12® is defined as Eq. (l52ll . Similarly, the right- 
hand-side of Eq. (IB21I) satisfies 


Therefore, Eq. (|B21ll is reduced to Eq. (|5T]l . The per¬ 
turbative solution of Eq. m with respect to the small 
inelasticity is given by Eq. dZQj). 

Similarly, we derive the differential equation for the 
thermal conductivity k. Multiplying Eq. (IB6|) by 
S{Vi)/T and integrating over Vi, we obtain 

^ (3c(°)«r) -h ly (/(°\^) 

= ^JdV,S{V,)-A{V,). (B24) 

The second term on the left-hand-side of Eq. (IB24I) is 
written as 

i Jdv,s(y,)j^^'> (/(°),^) = (B25) 

where 12® is given by Eq. (1551) . The right-hand-side on 
Eq. (IB24P satisfies 

If 15 nT 

- / dV,S{V^) ■ A{V^) = (1 + 2a2). (B26) 

1 J 2 m 

It should be noted that terms proportional to o„ (n > 3) 
vanish due to the orthogonality of the Sonine polynomi¬ 
als. Therefore, Eq. (IB24I) is reduced to Eq. (l53l) . The 
solution of Eq. (1551) is given by Eq. (171)) . 

Similarly, multiplying Eq. (IB7|) by S{Vi)/T and inte¬ 
grating over Vi, the coefficient p, is given by Eq. (ESI). 


J dViD,^{V)Cij{Vi) = lOnT. (B23) 

I 

Appendix C: Determination of a2 and as 


In this appendix, we determine the coefficients 02 and 03 using the moments of the dimensionless collision integrals 
It is useful to introduce the basic integral 


Jk,l,m,n,p,c.= / dC / dCi 2 / dfed(x, C12, 0 I ^12 ' fc |'+“<^(C')</.(Ci2 )C'=c‘ 3 (C • 012 )"* (C • fc)" (Ci2 • fcf, (Cl) 


with C = (ci -I- C2)/2. This is rewritten as 

dk^l,m^n,p^a — 2 


/ fc + m + n + 3 \ ^ /n\ ^ Tdesin^+i© 

V 2 / ,=n VJ/ r(-^)7o 


X / dc 


12 


db be 


j=o 

l+m+p+a +3 21 


Sin ■ 


X 


cos^ - exp ( “2^12 


I cos™+"-2 0 

(C2) 
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For a = 0 and n = 0, 1 and 2, Eq. (IC2I) reduces to 


o — 3)/2 

=-—j- [1 + r 


Jk,l,m, 2 ,p ,0 — 


m + 1 

2 ^ — {k-\-m—2)/2 


m + 2 

2”(fc+m —1)/2 


[i-(-ir]r 


fc + TO + 3 


fc + TO + 4 


dci 2 

dci 2 


db be 


l+m+p +3 . p X 


12 


Sin'- - exp --Ci2 , 


1 


db be 


l+m+p +3 „• p+1 X 


sin^ -exp --Ci2 


[i + (-i)™]r 


A: + TO + 5 


(to + 1)(to + 3) 

X J dci2 J sin^ ^ ^1 + TOsin^ exp ^ 


2^12 I > 


(C3) 

(C4) 


(C5) 


respectively. These integrals recover the previous results in the hard core limit [^ . In this paper, we only consider 
the nearly elastic case 1 — e 1. We assume that the coefficients 02 and 03 are proportional to 1 — e. When we use 
the truncated distribution function Eq. (l56l) . we rewrite the n-th moment Mp = — f dcicf {p G N) as 


Mp = 


dCdci2dkd{x,Ci2,0Wi2 ' k\(l>{ci)4>{c2){ci2 ■ kf 

X [1 + a 2 {S 2 (cf) + S 2 {c 2 )) + a3(S3{Ci) + S' 3 (c 2 ))] A (c^ + c^), 


(C 6 ) 


where, we have ignored the terms proportional to af, a§, and 0203 , because they are the order of (1 — e)^. The explicit 
forms of A(c^ + c^) for p = 2, 4, and 6 are, respectively, given by 


A(ci + C2) = - e0(i>max - by'' 


cos^ 0 
cos^ 6 


(ci2-fc)2+0(e2). 


(C7) 


A(ct + 4)=- 8(C ■ ci2)(C • fe)(ci2 • k) + 8(C ■ fc)2(ci2 • k^ 


-j- 6 0(&max 


,cos^ 0 
cos^ 9 
1 


- 2 C^(ci 2 • ky - -cf2(ci2 • ky + 4(C • C 12 )(C • fc)(ci2 • k) - 8{C ■ kyici 2 ■ ky 

A(c® + cl)=- 24C2(C • C 12 )(C • fc)(ci 2 • k) + 2AC^{C ■ fe)2(ci2 • k)^ 

- 6 c? 2 (C • C12)(C • fc)(ci2 • fc) + 6 c? 2 (C • fc) 2 (ci 2 • fc)" 

, cos^ 0 


+ 0(e"), (C 8 ) 


+ e0(6max - by-^ 


cos^ 9 


3 C 4 (ci 2 • fc)" + 2 <^"c?2(ci 2 • fc)" - 12 C"(C • Ci2)(C • fc)(ci2 • fc) + 24 C'"(C' • fc)"(ci2 • fc)" 

+ 44(C12 • ky - 3c?2(C • C12)(C • fe)(Ci2 • fc) + 6c?2(C' • fe)"(ci2 • fc)" + 3 (C • C12)"(C12 • ky 
lb 


-12{C ■ C 12 )(C • fc)(ci2 • ky + 12{C ■ fc)"(ci2 • ky + 0(e"). 


(C9) 


Then, we explicitly write A 42 , M 4 , and as 


M2 = \f2y (S'! + 02^2 + OSS's) , 
A44 = \/27r (Ti + 0272 + a^Ty , 
Mq = (Z?i + 02772 + 03773 ), 


(CIO) 
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where 


S'! =e / dci2 I dbb{iy -b )cl2e^p[-7:012] +0{e ), 


1 


0 ^0 

00 nbrr 


1 


S2=e—J dci2 j - 5 )ci 2(15 - 10ci2 + Ci2)exp |^--Ci2j +0(e ), 

1 pCO pb^nax ^ ^ ^ X T 

*^3 =£ 77 :;: / c;ci 2 / dhb{v'^ - 6^)ci 2(105 - 105ci2 + 21 ci 2 - Ci 2 )exp ( --c ?2 ) + C>(e^), 


Ti =e- 


192 JO 
1 

Jo 


dci2 / dhb{v^ - 6 ^)c? 2(5 + cf2)exp --0^2 +C»(e ), 


72 =^^ c^^^Ci2sin^X^°^ exp 


1 rOO POjxiax ^ ^ ^ X 1 

— J dci2 J - 6^)ci2(-25 - 23ci2 - 5ci 2+ Ci2)exp |^--Ci2 

/*^max __ _ __ / 1 

+ / 7 ci 2 / db b{v^ — b‘^)c [2 sin^ — 7 — exp I —-C 12 


+ i / (ici 2 db 6 c[ 2 X^^^ sin^ exp 

73 y dci2 J 7&6ci2(7-Ci2)sin^X^°^exp ^-ici2^ 


+ 0(e"), 


^ ^00 pOxnax / 1 

— J dci 2 J dbb{iy‘^ - &^)ci2(-525 - 168ci2 - 6 ci 2 + 16ci2 - 0 * 2 ) exp (^- 2^12 

1 /*°^ r^raax _ _ / 1 

+ - / dci 2 dbb{iy'^ - b‘^)cl2{7 - cl 2 )sm‘^-—exp i--cj. 


+ ^y dci2 J 7&&c[ 2(7-c?2)x^^^siii2 2x^^Yxp 


+ 0(e"), 


77i =e 


PCO pbxaax ^ ^ ^ 1 \ 

J dci2 j - 6 ^)c^ 2(35 + 14 c? 2 +ct 2 )exp (^--c^2j+C’(e^), 

7I2 y 7 ci 2 y db 6 ci 2(7 + C12) sin^ x^°^ exp ^-^€12^ 


256 


7ci2 / dbb{v^ — &^)ci2 (—595 — 252 cj 2 — 180^2 + 4ci2 + C12) exp I —-c' 


-12 


O ^00 ^Omax _ _ ^ (Uj / 1 

+ 7/ dci2 rf 65 (i/^-& 2 )ci 2 ( 7 -Ci 2 )sin^—-exp(--ci2 


/o 


/o 


O ^00 /-hmax _ _ _ _(0) / 1 

-7/ dci2 - &2)c?2siii'‘-—exp (--c' 


/o 


/o 


12 


3 

16 


y fici2 y rf6642(7 + Ci2)x^^^sin2 2x^°Yxp^-ic?2^ 


+ 0(e"), 


(Cll) 

(C12) 

(C13) 

(C14) 


(C15) 


(C16) 

(C17) 


(C18) 



18 


Ds J dci2 db &Ci2(35 - 0^2) sin^ exp ^-^012^ 


1024 


16 


dci2 / db b{u^ — 6^)ci2 (—5145 — 1785 ci 2 + 798 ci 2 + 22ci2 + 7ci2 — exp ( — 


^0 

' P^max ^( 0 ) 

dci2 / d 6 &(i/^ — 5^)012(35 — 14 ci 2 + C12) sin^——exp ( 


:„2 

A 12 


o pOO pOmax _ _ _ / 1 

+ g y dci2 J - 6 ^)c?2 (7-Ci2)siii'‘— exp |^--Ci2 

+ ^ / rfci2 J db 6^2(35 - 0^2)%^^^ sin^ 2x^°^ exp 


0{e^). 


(C 19 ) 


Here, we only show the lowest order of each term. Here, AI4 and A 4 e are, respectively, related to A42, the fourth 
moment (c^) and the sixth moment (c®) as 


^M2 (c'^) = Mi 
2M2 (c®) = Me 


(C 20 ) 


Substituting Eqs. (IClOp into Eq. (IC 20 I 1 with (c"^) = ( 15 / 4)(1 + 02) and (c®) = ( 105 / 8)(1 + 3a2 — 03), we obtain the 
simultaneous equations with respect to 02 and 03 as 


' (55i + 552 - T 2 ) 02 + (553 - Ts) 03 = Ti - 55i 

315 105 \ 105 105 \ 105 

^-5i H —^‘^2 — 112 j 02 + 1- — D 3 j as = Di - 


(C 21 ) 


These equations can be solved easily and the explicit forms of 02 and 03 up to e order are given by Eqs. (EZD-dnH). 
Thus, we explicitly write A42, A44, and Me up to the first order of e as 


A 42 =e\f^ / dc 


12 


dbb{v^ - 6^)cj2exp ( “2^12 ) + 


Mi =€ 


Me =e 


,/Ott rbma,x _ / 1 , 

(l-e)y dci 2 J d6 6 (z /2 - 6^)c® 2 (5 + C12) exp f--C' 

J d 6 6 ci 2 sin 2 x^°^ exp ^-ici 2 ^ 

+a3^^-^y d6^2 (7-^ 2 ) sin^X^°^ exp 


+ 0 (e% 


16 


(1 — e) / dci2 / d6 5 (z/^ — 6 ^)c®2 (35 + 14 cj 2 + C12) exp I —-c' 


(i)3\/^ 

+“2 16 

(1) 3\/^ 


-12 


H”Ci 


64 


J dci2 J db 6ci2 (7 + c^a) sin^ exp ^-^012^ 
J dci2 J db 6ci 2 (35 - C12) sin^ x^”^ exp 


0{e^). 


(C 22 ) 


(C 23 ) 


(C 24 ) 
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Appendix D: Calculation of Q’L and fl® 


In this appendix, we calculate the quantities 11® and H® in Eqs. 


CiCj C Sij/3^ Dij{c2^^ -^^ 1 (^ 2 ) 


V 

is rewritten as 


and ([55]). From the definition, Dij{c) = 


Dij{c2)A D^j{cl) + bij{c2) 


— [ ^2i^2j Q 


V'^2 


^2i^2j C-iiCij C2iC2j ij {ci + C2 C2) 

1 


=cf2(C • fe)(ci2 • k) - -C^2(ci2 • kY - 2{C • Ci2)(C • fe)(ci2 ■ k) + {C ■ C^2){ci2 ' kY 
+ 2 {C ■ kY{ci 2 ■ kY - 2 {c ■ fc)(ci2 • kY + i(ci2 • kY 


i0(6max 


, COS' 


2 0 


cos^ 0 


^C'^(ci 2 ■ kY - ^Ci 2 (C' • fc)(ci 2 • fc) + ^Ci 2 (ci 2 ■ kY + (C ■ Ci 2 )(C • fe)(ci 2 • fc) 
-he ■ C 12 )(C 12 • kY - 2(C • kY{ci 2 ■ kY + 2{C ■ fc)(ci 2 • kY - i(ci 2 • kY 

D 2 


+ 0{eY. 

Substituting this result into Eq. (l5^ . we obtain Eq. (1651) . 
For n®, we rewrite ^( 02 ) • A S{ci) + S{c 2 ) as 


(Dl) 


S{c 2 )-A S{cY + S{c 2 ) 


= {^2 - [(c'l ■ C2) cf + (4 • C2) C2^ - (ci • C2) cl - (C2 • C2) C2] 


=C\lYC ■ fc)(Ci 2 • fc) - 4C2(C • Ci 2 )(C • fc)(Ci 2 • fc) + CYC ■ Ci 2 )(Ci 2 • fc)^ + • fc)2(ci2 • fc)^ 

- 2 C 2 (C • fe)(ci2 • kY + \c\YC ■ fc)(ci2 • fc) - 2 c 22 (C • Ci2)(C • fc)(ci2 • fc) + ^clYC ■ C12)(C12 • fc)" 

+ clYC ■ fe)"(ci 2 • kY - \clYC ■ fc)(ci 2 • kY - ^clYC ■ fc)(ci 2 • fc) + 4(C • C 12 )"(C • fc)(ci 2 • fc) 

- (C • C12)"(C12 • kY - 4 (C • C12)(C • kY{ci 2 ■ kY + 2{C ■ Ci2)(C • fc)(ci2 • kY + 10 (C • C12)(C • fc)(ci2 • fc) 

- ^{C ■ C 12 )(C 12 • fc)" - 10(C • fc)"(ci 2 • fc)" + 5(C • fc)(ci 2 • kY 


+ e0(&max - by- 


cos 


20 


cos^ 0 


-CHc ,2 ■ kY - Ic^Uc ■ fe)(ci 2 • fc) - lcyiYc ^2 ■ kY + 2CYC ■ c, 2 )iC ■ fe)(ci 2 • fc) 


+ CYC ■ C 12 )(C 12 • kY - 4C2(C • fc)2(ci2 • fc)" + 2C2(C • fc)(ci 2 • fc) + -C2 (ci 2 • fc)" 

-■^cf 2 (C' • fc)(ci 2 • fc) + c 22(C • Ci 2 )(C' • fc)(ci 2 • fc) - c22(C • fe)"(ci 2 • fc)" + ^ 012(0 • fc)(ci2 • fe)^ 

+ ^c 22(C • fc)(ci 2 • fc) - 2(C • C 12 )"(C • fe)(ci 2 • fc) + 4(C • C 12 )(C • kY {^12 ' ^Y 

-2{C ■ C 12 )(C • fe)(ci 2 • kY - 5{C ■ C 12 )(C ■ fe)(ci 2 • fc) + 10(C • fc)"(ci 2 • fc)" - 5(C • fe)(ci 2 • ky 


+ 0(e"). 


(D2) 


Substituting this into Eq. dSSj), we obtain Eq. (Ifi5ll after the long and tedious calculation. 


Appendix E: High and low temperature expansions 

We can evaluate the explicit forms of the transport co¬ 
efficients in terms of high temperature expansion. We can 


also evaluate the dissipation rate AI 2 as a low tempera¬ 
ture expansion, though it diverges in the low temperature 
limit. 

First, we discuss the high temperature expansion. 
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From Eq. m, we expand v as 


Similarly, k and fj, are, respectively, expanded as 


ly = Jl + 


2 e 

W 2 




for T/e 1. Substituting Eq. (IE1|) into Eq. (IS^ . we 
expand A ^2 in terms of the small parameter e/T as 


Mi°^ = 0, Mi^'’ = 




O 


if} 


(E2) 


with 

= 2^, = 2\/^. (E3) 

Similarly, and fl® are, respectively, expanded as 


^e(O) ^ ^e(O.O) ^ £.f^e(0,l) ^ ^ (^4) 

= + (E5) 

^e(0) ^ ^e(0.0) ^ £.f^e(0,l) ^ ((|;)') , (E6) 

=f2f’°)+o(|;) (E7) 

with 

Of ■«) = -4v^, Of’0) = (E8) 


of-i) = ^(A - 1) {2(15A^ + 15A3 + 2A2 + 2A + 2) 
+3A^(A + 1)(5A^ — 1) log , (E9) 


^e(0,0) ^ f|e(1.0) ^ _ 


1989v^ 
320 ’ 


(ElO) 


oe(o.i) ^ _ 1 ) {2(15A^ + 15A3 + 2A2 + 2A + 2) 


24 


+3A2(A + l)(5A2-l)log 


A- 1 
A + 1 


• (Ell) 


Next, let us calculate the expansions of the transport 
coefficients. Substituting Eqs. (IE2])-(|E5]) into Eqs. dZlD 
and (1731) . we expand rj as 


+ |T»“’'‘> + of)"), (E12) 

,(i> = ,(i-»)+o(i) (E13) 

with 


^(0.0) 

^( 0 . 1 ) 


16d2 

= v' 


TT ’ ' 3840 ' 


(E14) 


= jo,o)^ —{2(15A'^ + 15A^ + 2A^ + 2A + 2) 


96 

+3A^(A + 1)(5A^- l)log 


A- 1 
A + 1 


. (E15) 



(E16) 

„(1) , + 

(E17) 

f.‘”> = 0. + 0 (£) 

(E18) 

with 


,+0.0) _ 75 rT~ _ 539 ^^(0 0) 

64c?2 V TTTO ’ 1280 ’ 

(E19) 

,,(0.1) ^ „(0.0) A - 1 ^ ^ 2A2 

96 

+ 2A + 2) 


+3A2(A + l)(5A2-l)log^|, (E20) 


n 1185 pT^ 

^ ’ “ 1024nd2 V 

Let us also calculate the low temperature expansion of 
Ad 2 - From Eq. ([6]), we expand 12 as 



V2ci2 T 
4 e 


+ 0 



(E22) 


Substituting Eq. (IE22I) into Eq. (IS^ . we can expand Ad 2 
in terms of the small parameter T je as 


M‘? = 0. + M'iS' + O (f) 

(E23) 

with 

Adf= 2V^A^ Adf°^ = 2\/^. (E24) 


Appendix F: DSMC algorithm 

In this appendix, we briefly summarize the DSMC pro¬ 
cedure |2J-[^, 1^ , which is a numerical technique to 
obtain the solution of the Boltzmann equation at t + At 
from that at t. For small At, the velocity distribution 
function at t + At is given by 

/(n,t +At) =/(t;,t) + ^^^L^At. (FI) 

Substituting the Boltzmann equation (13 into this, we 
obtain 

f{v, t + At) = (1 — AtD + At J) f{v, t) 

= (1 + At J) (1 - AtD) f{v, t) + 0 (Af ), 

(F2) 

where we have introduced Df = n • V/ and Jf = I{f, /) 
for simplicity. Equation (IF2I) shows that the time evo¬ 
lution of the velocity distribution function can be sep¬ 
arated into two parts: advective process and collision 
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process. According to this separation, DSMC iteration 
is as follows: (i) We determine the time step At smaller 
than L/vmax, where L is the system size and Wmax is 
the maximum speed among the particles, which is eval¬ 
uated as Umax = 5 ut with the thermal velocity ut- In 
this paper, we adopt At = 0.2L/umax. (h) We move 
the particles during At without any collisions. This cor¬ 
responds to update the distribution function f*(v,t) = 
(1 — AtD) f{v,t). (hi) We modify the velocities of the 
particles due to collisions. We randomly determine the 
collisions without taking into account the actual positions 
of the particles. The square of the collision parameter, b^, 
is chosen in the range 0 < 6^ < X^cP at random. A pair of 
colliding particles change the velocities according to rule 
in Eqs. ra and (|2n)l for a hard core collisions and Eqs. 
(HU) and ( 1211 ) for a grazing collision Here, the number 
of collisions Nc is evaluated as 7r(Ad)^iV^UniaxAt, which 
is proportional to the total cross section, the maximum 
speed, and the time step At. This process corresponds to 
obtain f{v, t + At) = (1 -I- At J) f*{v, t). (iv) We update 
the time t + At. 


Appendix G: Estimation of the trapping 
temperatnre 

As stated in Sec. m we ignore the trapping process 
by the potential well throughout the paper. In this Ap¬ 
pendix, we briefly discuss the critical condition which val¬ 


idates this approximation using a simple one-dimensional 
model. Let us consider a process in which two particles 
approach from far away relative speed v in the frame 
that the target is stationary. When the particle enters 
the potential region {d < r < Ad), the velocity becomes 

Vin = Jv'^ + — (Gl) 

V m 

from the energy conservation {l/2)m,.v^ = {l/2)mrvf^—e 
with the reduced mass mr = m/2. After the inelastic 
scattering on the hard core {r = d), the velocity changes 
from Vin to —euin. When the particle is trapped by the 
potential, the energy is negative, that is, (l/2)mre^v/^ — 
s < 0. Using Eq. (IGII) . the trapping condition is given 
by 

^^<^^trap = \/-(l-e)+0(l-e). (G2) 

V m 

The corresponding granular temperature is given by 

T < Ttrap = imi^Lp = 4e(l - e) + O ((1 - e)^) . (G3) 

The critical trapping temperature is given by Ttrap — 
O.Ode for e = 0.99, which is consistent with the lower 
bound of our theory as shown in Figs. Eland El 
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